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20.  ABSTRACT  (Continued) 


flutter  and  divergence  flight  envelope  and  may  lead  to  fatigue  damage  of  the 
system  even  though  aeroelastic  instability  is  not  encountered. 

During  a  previous  AFOSR  sponsored  study,  an  asymptotic  expansion  approxi¬ 
mation  technique  was  developed  to  predict  the  limit  cycle  response  of  aerodynamic 
surfaces  with  discrete  structural  nonlinearities.  As  a  part  of  this  study, 
numerical  simulation  results  were  used  to  evaluate  the  adequacy  of  the  asymptotic 
expansion  solutions  to  the  nonlinear  problem.  A  "simple"  trapezoidal  numerical 
integration  procedure  was  used  to  obtain  these  simulation  results.  Results  of 
this  previous  study  uncovered  shortcomings  with  the  trapezoidal  integration 
scheme  Based  on  these  results  it  was  concluded  that  the  application  of  more 
refined  integration  techniques  to  the  nonlinear  aerodynamic  surface  problem 
needed  to  be  investigated. 

Thus  the  present  study  was  undertaken  with  the  objective  to  evaluate,  on  a 
comparative  basis,  different  numerical  simulation  approaches  for  predicting 
limit  cycle  response  of  aerodynamic  surfaces  containing  discrete  structural 
nonlinearities.  Results  from  such  simulations  are  needed  to  compare  and  evaluate 
approximate  solutions  for  the  limit  cycle  response  of  nonlinear  systems.  In 
addition,  these  simulation  results  provide  information  concerning  the  nature  of 
the  nonlinear  system  response  which  may  be  used  to  aid  in  understanding  the 
mechanism  of  the  aerodynamic  surface  dynamics  and  in  understanding  the  response 
of  nonlinear  systems  in  general. 

* 

Three  numerical  integration  techniques  were  selected  for  evaluations  (1) 
fourth-order  Runge-Kutta,  (2)  eighth-order  Shanks,  and  f-3)  fourth-order  Adams- 
Moulton  predictor-corrector.  The  results  of  the  three  simulation  techniques 
compared  well  with  each  other.  In  addition  they  yielded  improved  correlation 
with  previously  developed  approximate  solutions  when  compared  to  the  correlation 
of  the  approximate  solutions  using  the  trapezoidal  integration  scheme.  It  was 
concluded  that  any  one  of  the  three  numerical  techniques  is  appropriate  for  use 
in  determining  limit  cycle  response  of  an  aerodynamic  surface  containing 
discrete  structural  nonlinearities.  From  a  computational  efficiency  point  of 
view',  the  Runge-Kutta  approach  appears  most  attractive  for  this  type  of  problem. 

In  spite  of  the  improved  results  obtained  with  the  numerical  integration 
techniques  evaluated  during  this  study,  there  remain  regions  where  the  correla¬ 
tion  between  numerical  and  approximate  prediction  is  inconsistent.  The  need  for 
additional  research  to  address  these  inconsistencies  has  been  identified.  This 
research  will  provide  an  improved  understanding  of  the  nonlinear  response  char¬ 
acteristics  of  an  aerodynamic  surface  containing  discrete  structural  nonlinear¬ 
ities. 


ACKNOWLEDGMENT 


This  report  was  prepared  by  the  Structural  Dynamics  and  Loads  Department, 
McDonnell  Douglas  Astronautics  Company-St.  Louis,  Missouri,  under  Contract 
F49620-84-C-0123,  "Investigation  of  an  Asymptotic  Expansion  Technique  to 
Analyze  Limit  Cycle  Response  of  Aerodynamic  Surfaces  with  Structural  Non- 
linearities."  The  program  was  administered  under  the  direction  of  the  Air  Force 
Office  of  Scientific  Research  with  Dr.  Anthony  K.  Amos  as  Program  Manager. 


*1? 


"s  •  >  i-  - 


<e.  ; 


■  '•f»J  la 
-12. 


3 


TABLE  OF  CONTENTS 


Page 


1.0  INTRODUCTION .  1 

2.0  NUMERICAL  INTEGRATION  TECHNIQUES .  10 

2.1  FOURTH-ORDER  RUNGE-KUTTA .  12 

2.2  EIGHTH-ORDER  SHANKS .  13 

2.3  FOURTH-ORDER  ADAMS-MOULTON .  15 

3.0  NUMERICAL  SIMULATION  RESULTS .  17 

3.1  RIGID  AERODYNAMIC  SURFACE .  17 

3.2  FLEXIBLE  AERODYNAMIC  SURFACE .  20 

4.0  CONCLUSIONS .  31 

5.0  REFERENCES .  35 

APPENDICES 

A.  AERODYNAMIC  SURFACE  CONFIGURATION .  37 

B.  SIMULATION  RESULTS  FOR  FLEXIBLE  AERODYNAMIC  SURFACE .  41 


v 


LIST  OF  ILLUSTRATIONS 


Figure  Title  Pag 

1  Models  of  Aerodynamic  Surface  Aeroelastic  Response .  2 

2  Type  of  Structural  Non]  iije  ar  it  Fes .  3 

3  Aerodynamic  Surface  Configuration .  5 


4  Root  Roll  Response  Using  Trapezoidal  Integration  Technique _  7 

5  Root  Pitch  Response  Using  Trapezoidal  Integration  Technique...  8 


6  Simulation  Results  for  Rigid  Aerodynamic  Surface  with  Root 

Roll  Preload  Nonlinearity .  18 

7  Simulation  Results  for  Rigid  Aerodynamic  Surface  with  Root 

Pitch  Preload  Nonlinearity .  19 


8  Root  Roll  Simulation  Results  for  a  Flexible  Aerodynamic 

Surface  with  Two  Preload  Nonlinearities . 

9  Root  Pitch  Simulation  Results  for  a  Flexible  Aerodynamic 

Surface  with  Two  Preload  Nonlinearities . 

10  Aerodynamic  Surface  Response  for  a  Dynamic  Pressure  of  50  psi 

11  Root  Roll  PSD  for  a  Dynamic  Pressure  of  50  psi . 

12  Frequency  Content  of  Root  Roll  Response  Versus  Dynamic 

Pressure . 

13  Lowest  Two  Frequency  Components  as  a  Function  of  Dynamic 

Pressure . 

14  Aerodynamic  Surface  Response  for  a  Dynamic  Pressure  of  70  psi 

15  Root  Roll  PSD  for  a  Dynamic  Pressure  of  70  psi . 

16  Comparison  of  Computational  Efficiency  for  Numerical 

Integration  Techniques . 

17  Region  of  Correlation  Anomaly . 

A1  Aerodynamic  Surface  Geometry . 

A2  Aerodynamic  Surface  Inertia  Properties . 

A3  Aerodynamic  Surface  Cantilerer  Modes . 


LIST  OF  SYMBOLS 


A 

B 

D 

h 

I 

K  ( X ) 

M 

mn 

P 

PF 

Q 

q 

qn 

s 

t 

x 

Y 

4> 

0 

W(|),U)Q 

U)n 


Amplitude  of  motion  at  surface  root  support 

Elements  of  aerodynamic  forcing  function  matrix 

State  variable  dynamic  matrix 

Numerical  integration  step  size 

Rigid  aerodynamic  surface  inertia 

Spring  rate  of  root  springs 

System  mass 

Flexible  aerodynamic  surface  generalized  masses 
Preload  parameter 

Aerodynamic  surface  participation  factors 
Aerodynamic  matrix  for  state  variable  formulation 
Dynamic  pressure 

Flexible  aerodynamic  surface  generalized  coordinates 
Deadspace 
T  ime 

Matrix  of  aerodynamic  surface  degrees  of  freedom 
Matrix  of  aerodynamic  surface  state  variables 
Displacement  in  root  pitch 
Displacement  in  root  roll 

Uncoupled  natural  frequencies  of  linearized  system 
Flexible  aerodynamic  surface  natural  frequencies 


Subscripts 
A. .  .K 

n 

<t> 

9 


Indicies  associated  with  intermediate  numerical  integration 
time  steps 

Indicies  associated  with  numerical  integration  time  step 
Associated  with  root  pitch  degree  of  freedom 
Associated  with  root  roll  degree  of  freedom 


1.0  INTRODUCTION 


Defining  the  flutter  and  divergence  characteristics  of  aerodynamic  sur¬ 
faces  is  a  basic  requirement  in  assuring  structural  and  performance  integrity 
of  a  given  design  for  its  operational  environment.  The  response  characteris¬ 
tics  of  the  divergence  and  flutter  phenomena  are  illustrated  in  Figures  1(a) 
and  (b).  For  systems  containing  structural  nonlinearities,  another  mode  of 
aeroelastic  response--limit  cycle  osci llation--may  be  present.  Limit  cycle 
response.  Figure  1(c),  is  characterized  by  steady  state  oscillation  whereas 
divergence  and  flutter  are  unstable  motions  with  increasing  amplitude.  The 
potential  of  limit  cycle  response  is  of  importance  since  these  oscillations  may 
occur  within  the  aerosurface  flutter  and  divergence  flight  envelope.  The 
amplitude,  frequency,  and  duration  of  these  potential  limit  cycle  oscillations 
are  of  interest  for  a  fatigue  assessment  of  the  system. 

Frequently  aerodynamic  surface  hardware  designs  do  have  nonlinearities  in 
the  surface  itself,  support  structure,  and/or  control  actuators  as  a  result  of 
manufacturing  tolerances,  design  characteristics,  and/or  freeplay.  When  these 
nonlinearities  exist,  the  assumption  of  a  linear  force-displacement  relation¬ 
ship  is  no  longer  justified.  In  this  situation  an  understanding  of  the 
nonlinear  effect  on  the  dynamic  behavior  is  required  to  evaluate  system 
response. 

The  effects  of  structural  nonlinearities  on  aeroelastic  analysis  have  been 
studied  both  analytically  and  experimentally.  Reference  1  through  9.  In  these 
studies  several  nonlinearities  that  are  typically  encountered  in  aerodynamic 
surface  designs  were  considered.  The  primary  thrust  of  the  early  work 
attempted  to  establish  a  basic  understanding  of  the  nonlinear  system  response 
employing  analog  computers.  The  later  analytical  studies.  References  4  and  6, 
employed  the  describing  function  techniques  to  characterize  the  nonlinear 
behavior  of  an  aerodynamic  surface.  A  majority  of  the  work  addressed  freeplay 
nonlinearities  of  the  type  shown  in  Figure  2.  These  nonlinearities  are 
representative  of  a  deadband  or  "slop"  in  the  root  support  structure  of  an 
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TYPES  OF  STRUCTURAL  NONLINEARITIES 


aerodynamic  surface  with  and  without  e  preload.  The  approach  discussed  in 
Reference  5  is  a  procedure  in  the  frequency  domain  which  iteratively  obtains  a 
consistent  solution  in  terms  of  response  amplitude  and  effective  stiffness  of 
the  nonlinear  structural  element. 

The  describing  function  approach,  References  4  and  6,  uses  a  one  term 
Fourier  Series  expansion  of  the  force  to  account  for  the  effect  of  the 
nonlinear  stiffness  on  aerodynamic  surface  response.  This  method  yields 
satisfactory  results  when  the  amplitude  of  displacement  is  greater  than  the 
magnitude  of  the  freeplay.  Figure  2(a),  or  freeplay  plus  preload.  Figure  2(b). 
However,  it  was  pointed  out  in  References  4  and  6  that  when  the  amplitude  of 
displacement  is  approximately  equal  to  the  magnitude  of  the  nonlinearity, 
significant  error  can  occur  as  a  result  of  neglecting  the  higher  harmonics  in 
the  series  expansion  of  the  force-displacement  relationship. 

An  analytical  study  sponsored  by  AFOSR  was  undertaken.  Reference  9,  to 
develop  an  improved  technique  for  predicting  limit  cycle  response  of  aero¬ 
dynamic  surfaces  with  discrete  structural  nonlinearities.  This  improved 
technique  was  to  retain  the  flexibility  of  the  describing  function  approach 
while  providing  greater  accuracy  and  generality  in  modeling  the  nonlinear 
system  behavior.  An  asymptotic  expansion  method  was  developed  to  model  the 
nonlinear  force-displacement  relationship  that  results  when  nonlinearities  of 
the  type  shown  in  Figure  2  are  introduced  at  the  aerodynamic  surface  support. 
The  primary  difference  between  the  asymptotic  method  and  the  describing 
function  method  is  the  capability  of  the  asymptotic  method  to  include  higher 
harmonics  in  the  representation  of  the  system  nonlinearity.  In  this  manner  one 
may  obtain  successively  higher  c.’der  approximations  to  the  limit  cycle 
response. 


Specifically,  the  problem  investigated  during  the  Reference  9  study  was 
the  limit  cycle  response  of  an  aerodynamic  surface  in  a  subsonic  airstream, 
Figure  3.  The  nonlinearities  shown  in  Figure  2  were  assumed  to  act  at  the  root 
support  springs  K<j>  and  Kg  of  the  structure.  This  problem  is  representative  of 
a  control  surface  with  a  loose  hinge  and/or  joint  slippage  in  the  surface 


support  structure  and/or  control  actuator.  The  aerodynamic  forces  acting  on 
the  surface  were  modeled  using  steady  state  aerodynamic  theory.  This  theory 
assumes  the  lifting  force  is  proportional  to  and  in  phase  with  the  torsional 
motion  of  the  surface  which  is  assumed  to  be  sinusoidal.  Simple  aerodynamics 
were  assumed  since  the  primary  objective  of  this  previous  study  was  the 
investigation  of  the  influence  of  structural  nonlinearities  on  aerodynamic 
surface  response.  As  discussed  in  Reference  6,  use  of  a  more  sophisticated 
aerodynamic  theory  can  substantially  change  the  linear  flutter  results 
employed  to  predict  nonlinear  aeroelastic  response,  but  has  no  impact  on 
interpretation  of  the  system  response  behavior  due  to  the  presence  of 
structural  nonl inearities. 

During  the  Reference  9  study,  numerical  simulation  results  were  used  to 
evaluate  the  adequacy  of  the  asymptotic  expansion  technique.  The  "exact" 
solutions  for  aerodynamic  surface  response  were  obtained  via  numerical 
integration  of  the  system  nonlinear  equations  of  motion.  These  "exact" 
solutions  were  then  compared  with  the  asymptotic  expansion  predictions  to 
assess  the  accuracy  of  these  predictions.  The  numerical  simulation  approach 
employed  during  this  previous  investigation  was  based  on  a  "simple"  trape¬ 
zoidal  integration  method.  Reference  10. 

It  became  apparent  during  the  Reference  9  investigations  that  there  were 
situations  when  the  trapezoidal  simulation  technique  exhibited  numerical 
stability  problems.  Examples  of  the  results  observed  during  the  Reference  9 
study  are  shown  in  Figures  4  and  5.  Here  simulation  results  employing  the 
trapezoidal  integration  technique  are  compared  with  the  asymptotic  expansion 
results  of  Reference  9.  This  example  is  for  a  flexible  aerodynamic  surface 
having  preload  nonlinearities  in  both  root  degrees  of  freedom.  The  differences 
in  simulation  results  when  compared  to  the  approximate  solution  may  be  noted. 
These  differences  indicate  that  the  trend  of  simulation  results  do  not 
correspond  to  those  of  the  asymptotic  expansion  technique.  Due  to  the 
"scatter"  of  the  simulation  results,  the  accuracy  of  these  results  was  in 
question. 
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Thus,  the  present  study  was  undertaken  to  evaluate  these  simulation 
shortcomings.  The  objective  of  this  study  was  to  evaluate,  on  a  comparative 
basis,  more  refined  numerical  simulation  techniques  for  predicting  the 
response  of  aerodynamic  surfaces  containing  discrete  structural  nonlinear¬ 
ities.  During  this  study  the  suitability  of  three  different  type  numerical 
simulation  techniques  were  compared  in  terms  of  predicting  the  large  amplitude 
limit  cycles  oscillations  of  a  nonlinear  aerodynamic  surface.  The  three 
numerical  techniques  investigated  were: 

•  Runge-Kutta 

.  Shanks 

.  Adams-Moul ton 

Details  of  these  numerical  integration  techniques  are  presented  in  Section 
2.0.  This  is  followed,  Section  3.0,  by  the  numerical  results  obtained  for  the 
three  numerical  procedures.  These  results  are  compared  with  both  the  data 
obtained  via  the  trapezoidal  integration  routine  and  the  asymptotic  expansion 
approach  of  Reference  9.  Study  conclusions  are  presented  in  Section  4.0. 


2.0  NUMERICAL  INTEGRATION  TECHNIQUES 


The  objective  of  this  present  study  was  to  evaluate,  on  a  comparative 
basis,  different  numerical  simulation  approaches  for  predicting  limit  cycle 
response  of  aerodynamic  surfaces  containing  discrete  structural  nonlinear¬ 
ities.  The  nonlinear  equations  of  motion  of  interest  are  of  the  form 

M  X  +  K ( X )  X  =  q  B  X  (1) 


These  system  of  equations  govern  the  nonlinear  aeroelastic  response  of  the 
aerodynamic  surface  conf igurat ion  shown  in  Figure  1.  The  detailed  elements  of 
Equation  (1)  are  given  as: 


For  this  study,  and  the  previous  investigations  of  References  4,  6  and  9, 
a  baseline  aerodynamic  surface  was  assumed  for  evaluating  the  coefficients  of 
Equation  (2).  This  baseline  configuration  was  based  on  the  Harpoon  Anti-Ship 
missile  quick-attach  control  surface.  Details  of  this  aerodynamic  surface  may 
be  found  in  Appendix  A. 

Three  numerical  integration  techniques  were  chosen  for  evaluation  during 
this  study.  These  were:  (1)  a  fourth-order  Runge-Kutta  method,  (2)  an 
eighth-order  Shanks  method,  and  (3)  a  fourth-order  Adams -Moulton  predictor- 
corrector  method.  Each  method  was  used  to  obtain  time  history  solutions  for 
the  equations  of  motion.  Equation  (1),  of  an  aerodynamic  surface  containing 
discrete  structural  nonl inearities. 
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All  three  methods  required  that  the  second  order  equations  of  Equation  (1) 
be  transformed  to  first  order  equations  using  an  appropriate  state  variable 
transformation.  Thus,  Equation  (i)  was  expressed  as: 

Y  +  DC V)  Y  =  Q  Y  (3) 


where 


.  \i 


The  matrices  D(Y)  and  Q  of  Equation  (4)  are  of  the  form 


D(Y )  = 


0  1  M_1K(X) 


:l  ;  0 


0  ;  q  M_1B 
0  '  0 


where  I  is  the  identity  matrix  and  the  remaining  terms  correspond  to  elements 
of  Equation  (1). 

The  following  notation  is  used  to  present  the  specific  form  of  each 
numerical  integration  technique.  For  the  initial  value  problem 

Y  -  f ( Y,  t),  Y ( t 0 )  =  Y0  (7) 

at  any  step  n  in  the  calculation,  the  available  quantities  are 

tn>  Y p  and  Y^  =  f(Yn,  t^)  (3) 


u 


Referring  to  Equation  (3)  the  quantity  f(Yn>  tn )  is  defined  as 

f(Yn,  tn)  =  Q  Yn  -  D(Yn)  Yn  (9) 

For  an  integration  step  size  set  equal  to  h  we  have 

tn+l  =  tn  +  h  (10) 

and  the  numerical  integration  procedure  computes  a  value  for  Yn+^.  The 
quanti ty 

Yn+1  =  f(Yn+1>  tn+i)  (11) 

is  then  calculated  and  the  integration  cycle  is  repeated. 

Each  of  the  three  numerical  integration  procedures  used  during  this  study 
are  described  in  detail  in  the  following  sections. 

2.1  Fourth-Order  Runge-Kutta 

The  Runge-Kutta  integration  method.  Reference  11,  was  selected  for  this 
study  because  it  is  a  basic,  widely  used  constant  step  integration  technique. 
This  method  uses  the  following  procedure  to  compute  Yn+^.  First,  the  sequence 
of  computations  given  below  are  performed: 

YA  --  Yn  +  (h /2)  Yn 

YA  =  f ( Ya,  tn  +  h/2) 

YB  =  Yn  +  (h/2)  YA  (12) 

YB  =  f  ( YB,  tn  +  h/2) 

YC  =  Yn  +  hYB 
Yc  =  f(YC’  ^  +  h) 
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Then  Yn+i  at  time  tn+i  is  defined  as: 


Yn+1  =  Yn  +  (h/6)  (Yn  +  2YA  +  2YB  +  Yc)  (13) 

2.2  Eighth-Order  Shanks 

The  eighth-order  Shanks  constant  step  integration  approach.  Reference  12, 
was  selected  as  a  means  to  incorporate  a  higher  order  method  of  the  Runge-Kutta 
type  approach  in  the  simulation.  This  was  of  interest  to  study  the  effects  of 
reducing  the  truncation  error  inherent  to  these  integration  methods.  To  obtain 
an  estimate  of  Yn+i,  first  the  following  sequence  of  calculations  are  made: 

YA  -  Yn  +  ( h/9)  Yn 

YA  =  HYA,  tn  +  h/9) 

YB  =  Yn  +  (h/24)  (Yn  +  3YA) 

yB  z  f(YB>  ln  f 

YC  =  Yn  +  (h/16)  (Yn  +  3YB) 

YC  =  f(YC,  tn  +  h/4)  (14) 

YD  =  Yn  +  ( h/500)  (29Yn  +  33YB  -  12YC) 

Yd  =  f ( Yq,  tn  +  h/10) 

YE  =  Yn  +  ( h/972 )  ( 33 Yn  +  4YC  +  125YD) 

Ye  =  f (YE,  tn  +  h/6) 
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Yf  =  Yn  +  (h/36)  ( -21Yn  +  76YC  +  125YD  -  162YE) 


VF  =  f ( Yp ,  tn  +  h/2) 

YG  =  Yn  +  (h/243)  (-30Yn  -  32YC  +  125Y0  +  99Yp) 

YG  =  <r( Yg,  tn  +  2h/3) 

Yh  =  Yn  +  (h/324)  ( 1175 Yn  -  3456YC  -  6250Yq 
+  8424Ye  +  242Yp  -  27Yg) 

YH  =  ^(YH,  tn  +  h/3) 

Y !  =  Yn  +  (h/324)  ( 293 Yn  -  852YC  -  1375YD  +  1836YE 

-  1 18Y  p  +  162 Yg  +  324 YH )  (14) 

CONTINUED 

Y i  =  f ( Yj ,  tn  +  5h/6) 

Yj  --  Yn  +  ( h/1620)  ( 1 303 Yn  -  4260YC  -  6875YD 
+  9990YE  +  1030YP  +  162Yj ) 

Yj  =  f(Yj,  tn  +  5h/6) 

YK  =  Yn  +  ( h/4428)  ( -8595 Y n  +  30720YC  +  48750YD 

-  66096YE  +  378Yp  -  729YG  -  1944YH  -  1296Y j  +  3240Yj) 
YK  =  Y(YK,  tn  +  h) 


The  value  for  Yn+i  at  a  time  tn+^  is  then  obtained  from: 


Yn+1  =  Yn  +  ( h/840)  (41Yn  +  216YE  +  272YF  (15) 

+  2  7  Yg  +  27 Yh  +  36Yj  +  180Yj  +  41 YK ) 

2.3  Fourth-Order  Adams-Moul ton 

The  Adams-Moulton  integration  technique.  Reference  11,  was  also  included 
in  this  study.  This  technique  was  selected  since  it  is  a  predictor-corrector 
procedure  and  thus  computationally  very  different  from  the  previous  two 
methods.  The  Adams-Moulton  procedure  uses  the  following  sequence  of  computa¬ 
tions  to  calculate  Yn+E.  The  predicted  value  at  tn+F  is 

pn+l  =  Yn  +  (h/24)  (55Yn  -  59Yn_j  +  37Yn_2  -  9Yn_3)  (16) 

The  corrected  value  at  tn+E  is  then  given  as 

Cn+1  =  Yn  +  (^24)  (9Pn+1  +  1 9 Yn  -  SY^  4  Yn_2)  (17) 

where 

pn+l  =  p(pn+l’  tn+i )  (18) 

and 

Yn+1  =  Cn+1  +  (19/270)  (Pn+i  -  Cn+1)  (19) 

The  predicted  and  corrected  values  are  used  to  evaluate  an  accuracy 
indicator  defined  as 

En  =  max[min(|Pn  -  Cp|,  |(P„  -  Cn)/C„| )] 


for  i  =  1,  2,  3,  .  .  ,  N 


(20) 


In  Equation  (20),  N  is  the  total  number  of  integrated  variables  and  i  denotes 
the  ith  integrated  variable. 

The  parameter  En  is  then  compared  to  two  constants,  Em-jn  and  Emax,  and  h 
i s  varied  as  fol lows : 

(1)  If  En_j  <  Em-jn  (j  =  0,  1,  2,  or  3);  h  is  doubled  and  the  integration 
is  restarted. 

(2)  If  Emi-n  <  En  <  Emax,  h  is  left  unchanged. 

(3)  If  En  >  Emax;  h  is  halved  and  the  integration  is  restarted. 

For  these  studies  Em-jn  was  set  at  5  x  10"9  and  Emax  at  5  x  10“6.  These  values 
were  selected  based  on  information  provided  in  Reference  13. 

As  can  be  seen  from  Equations  (16)  and  (17),  the  Adams-Moulton  method  is 
not  self-starting.  Each  time  the  integration  is  started  (restarted),  the 
Runge-Kutta  technique  is  used  to  compute  the  first  three  points.  The 
Adams-Moulton  procedure  is  then  used  until  a  restart  becomes  necessary. 
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3.0  NUMERICAL  SIMULATION  RESULTS 


Numerical  simulation  results  were  obtained  for  the  aerodynamic  surface 
nonlinear  equations  of  motion.  Equation  (1),  for  each  of  the  three  integration 
techniques  discussed  in  Section  2.0.  The  results  from  each  numerical  approach, 
for  a  uniform  set  of  physical  parameters,  were  compared  with  each  other  to 
evaluate  the  relative  accuracy  in  predicting  the  system  steady  state  limit 
cycle  response.  These  numerical  solutions  were  obtained  for  both  a  rigid  and  a 
flexible  aerodynamic  surface.  However,  a  majority  of  the  study  effort  was 
directed  toward  a  flexible  aerodynamic  surface  with  two  preload  nonlinear¬ 
ities.  This  configuration  was  of  primary  interest  since  the  results  of 
References  4  and  9  concluded  that  the  poorest  correlation  between  the 
approximate  solutions  and  simulation  results  was  observed  for  this  case. 

Power  Spectral  Density  (PSD)  analysis  studies  were  also  performed  to 
determine  the  frequency  content  of  the  time  history  simulation  results.  This 
frequency  content  information  is  of  interest  when  interpreting  the  numerical 
simulation  results.  As  discussed  in  Section  3.2,  review  of  the  frequency 
information  provided  insight  into  the  changing  nature  of  the  nonlinear  system 
response  as  a  function  of  varying  dynamic  pressure. 

3. 1  Rigid  Aerodynamic  Surface 

Results  of  the  numerical  simulations  for  the  rigid  aerodynamic  surface  are 
shown  in  Figure  6  and  7.  Shown  here  are  the  numerically  predicted  limit  cycle 
amplitudes  of  response  as  a  function  of  dynamic  pressure.  For  both  cases,  it 
was  assumed  that  the  uncoupled  root  roll  frequency  (ojq)  was  60  Hz  and  the  pitch 
frequency  )  was  215  Hz.  The  results  in  Figure  6  are  for  a  preload 
nonlinearity  in  the  root  roll  degree  of  freedom,  while  Figure  7  results  are  for 
a  root  pitch  preload  nonl ineari ty.  In  each  case  the  nonlinearity  was  defined 
by  a  deadspace  (S)  of  0.2  degrees  and  a  preload  (P)  of  0.1  degrees.  Thus  the 
S  over  P  ratio,  see  Figure  2,  was  two  for  both  cases. 


o  RIGID  AERODYNAMIC  SURFACE 

o  ROLL  PRELOAD  NONLINEARITY 

S q  =  0.2  deg 
P0  =  0.1  deg 

o  UNCOUPLED  FREQUENCIES 

oj.  =  215  Hz 
<P 

ujq  =  60  Hz 


DYNAMIC 

PRESSURE 

AMPLITUDE  RATIO 

(S/A) 

(psi ) 

RUNGE-KUTTA 

SHANKS 

ADAMS-MOULTON 

40 

0.34357 

0.34357 

0.34357 

60 

.31363 

.31363 

.31363 

30 

.31139 

.31139 

.31139 

90 

.11181 

.11181 

.11181 

95 

.05171 

.05165 

.05164 

o  RIGID  AERODYNAMIC  SURFACE 

o  PITCH  PRELOAD  NONLINEARITY 

S*  =  0-2  deg 
P^  =  0.1  deg 

o  UNCOUPLED  FREQUENCIES 

u).  =  215  Hz 
0 

oj  q  =  60  Hz 


DYNAMIC 

PRF^IIRF 

AMPLITUDE  RATIO  (S/A) 

i  1  \  L  J  J  Ul\  L 

(psi ) 

RUNGE-KUTTA 

SHANKS 

ADAMS-MOULTON 

40 

0.35331 

0.35295 

0.35295 

50 

.27660 

.27682 

.27682 

70 

.15124 

.15131 

.15243 

90 

.02935 

.02940 

.02944 

As  can  be  seen  by  the  data  presented  in  Figures  6  and  7,  the  results  from 
the  three  integration  procedures  compared  very  well  with  each  other  and  yield 
essentially  identical  results.  It  should  be  noted  that  the  computational  time 
was  quite  different  for  the  three  methods.  This  point  is  discussed  in  more 
detail  in  Section  4.0. 

3.2  Flexible  Aerodynamic  Surface 

The  numerical  simulation  results  for  a  flexible  aerodynamic  surface  are 
presented  in  Figures  8  and  9.  These  simulation  results,  amplitude  ratio  versus 
dynamic  pressure,  are  also  compared  to  asymptotic  expansion  results  of 
Reference  9.  The  Reference  9  simulation  results,  employing  the  trapazoidal 
procedure  are  also  shown  for  comparative  purposes. 

The  results  presented  in  these  figures  correspond  to  an  aerodynamic 
surface  with  preload  nonlinearities  in  both  root  degrees  of  freedom.  In  each 
case  the  deadspace  (S)  is  0.2  degrees  and  the  preload  (P)  is  0.1  degrees.  Thus 
the  preload  ratio  (S/P)  is  two  for  both  nonlinearities.  In  addition,  the 
uncoupled  root  roll  frequency  (coq)  was  60  Hz  and  the  pitch  frequency  (w^)  was 
215  Hz. 

Several  points  are  evident  from  the  data  shown  in  Figure  8  and  9.  First, 
the  three  integration  techniques  yield  similar  magnitudes  of  limit  cycle 
amplitude  ratios  for  a  given  set  of  physical  parameters.  Additionally,  this 
similarity  of  results  is  consistent  over  a  wide  range  of  dynamic  pressures  and 
associated  large  range  of  limit  cycle  response  amplitudes.  In  general  the 
results  obtained  with  the  refined  simulation  techniques  show  improved  cor¬ 
relation  with  the  asymptotic  expansion  results  when  compared  to  the  results 
obtained  via  the  trapezoidal  integration  scheme. 

However,  there  remain  regions  where  the  comparison  between  the  refined 
numerical  simulations  and  approximate  solutions  is  inconclusive.  To  more 
fully  evaluate  these  regions  of  inconsistent  results,  frequency  analyses  were 
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performed  on  the  simulation  developed  time  history  response  data.  This  was 
done  in  an  attempt  to  more  completely  define  the  mechanism  of  aerodynamic 
surface  dynamic  response  in  these  regions. 

These  investigations  of  response  frequency  content  were  done  by  first 
performing  PSD  analyses  of  the  root  degrees  of  freedom  displacement  data.  As  an 
example,  time  history  data  for  the  root  degrees  of  freedom  and  a  dynamic 
pressure  of  50  psi  are  shown  in  Figure  10.  The  corresponding  roll  response  PSD, 
(radians^/Hz  versus  Hz)  is  shown  in  Figure  11. 

From  data  such  as  shown  in  Figure  11,  the  frequency  content  of  the  system 
response  was  obtained.  This  was  accomplished  by  noting  the  frequencies 
associated  with  the  peak  or  highest  values  of  PSD.  A  summary  of  these  results 
is  given  in  Figure  12  for  response  in  the  root  roll  degree  of  freedom  over  a 
wide  range  of  dynamic  pressures.  The  highest  two  frequency  components  remain 
essentially  constant  over  the  range  of  dynamic  pressures  evaluated.  As  may  be 
seen  from  Figure  12,  there  is  significant  change  in  the  lowest  two  frequencies 
with  changing  dynamic  pressure.  Time  history  information  from  which  the  data 
presented  in  Figure  12  were  derived  are  presented  in  Appendix  B. 

These  lowest  two  frequencies  are  plotted  as  a  function  of  dynamic  pressure 
in  Figure  13.  Note  that  these  frequencies  are  the  average  values  for  the  root 
pitch  and  root  roll  response.  There  was  some  slight  difference  in  the 
frequencies  obtained  from  the  roll  and  pitch  PSD’s.  The  trend  shown  in  Figure 
13  is  similar  to  the  classic  frequency  coalescence  of  a  linear  system  flutter 
analysis. 

It  is  of  interest  to  note  that  there  is  a  change  in  the  nature  of  the  limit 
cycle  response  as  the  dynamic  pressure  approaches  the  linear  system  flutter 
boundary.  This  change  in  response  characteristics  may  be  seen  by  comparing  the 
results  shown  in  Figures  14  and  15,  for  a  dynamic  pressure  of  70  psi,  with  that 
of  Figures  10  and  11  for  a  dynamic  pressure  of  50  psi.  As  can  been  seen  in 
Figure  14,  the  root  displacement  has  developed  a  beat  type  characteristic. 
This  is  further  manifested  by  the  PSD  results  of  Figure  15.  From  the  data  in 


23 


o  FLEXIBLE  AERODYNAMIC  SURFACE 
o  TWO  PRELOAD  NONLINEARITIES 
s4>  ■  se  *  0-2  deg 

p4)  =  pe  =  °*1  de9 
o  UNCOUPLED  FREQUENCIES 

oi^  =  215  H2 
ujg  =  60  Hz 


TIME  (SEC) 


FIGURE  10  AERODYNAMIC  SURFACE  RESPONSE  FOR  A  DYNAMIC  PRESSURE  OF  50  PS  I 
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FIGURE  12  FREQUENCY  CONTENT  OF  ROOT  ROLL  RESPONSE  VERSUS  DYNAMIC  PRESSURE 
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FIGURE  15  ROOT  ROLL  PSD  FOR  A  DYNAMIC  PRESSURE  OF  70  PSI 


this  figure  it  appears  that  the  lower  two  frequencies  are  merging  toward  a 
single  value  within  the  accuracy  of  the  PSD  plot.  This  is  in  fact  the  result 
which  is  indicated  by  the  plot  in  Figure  13. 

Referring  to  Figure  9,  it  is  noted  that  the  comparison  between  the  three 
numerical  simulation  techniques  and  the  approximate  solution  is  inconclusive 
in  the  dynamic  pressure  range  of  55  to  65  psi.  In  fact,  the  numerical  results 
indicate  a  constant  amplitude  ratio,  or  "flat  spot",  over  this  range  of  dynamic 
pressures.  This  region  of  dynamic  pressures  precedes  the  changing  of  the 
response  characteristics  to  that  of  a  beat  phenomenon.  The  implication  of  this 
response  characteristic  is  discussed  in  Section  4.0. 
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4.0  CONCLUSIONS 


The  objective  of  this  study  was  to  evaluate,  on  a  comparative  basis, 
several  different  type  numerical  simulation  techniques  for  predicting  the 
limit  cycle  response  of  an  aerodynamic  surface  with  structural  nonlinearities. 
This  objective  was  of  interest  since  earlier  studies.  Reference  9,  uncovered 
shortcomings  with  a  "simple"  trapezoidal  integration  scheme.  Numerical  data 
obtained  with  this  trapezoidal  approach  were  used  for  comparison  with  and 
evaluation  of  approximate  solutions  to  the  nonlinear  problem.  Based  on  these 
previous  results,  it  was  felt  that  the  application  of  more  refined  integration 
techniques  for  nonlinear  aerodynamic  surfaces  needed  to  be  evaluated.  Thus  the 
present  study  was  undertaken. 

The  objective  of  this  study  has  been  met.  Limit  cycle  response  predictions 
were  obtained  for  the  aerodynamic  surface  including  discrete  structural 
nonlinearities  shown  in  Figure  1.  These  predictions  were  made  employing  the 
following  three,  quite  different  numerical  integration  schemes: 

•  Runge-Kutta 

.  Shanks 

.  Adams-Moulton 

The  results  of  the  three  simulation  techniques  compared  well  with  each  other. 
Thus  it  is  concluded  that  any  one  of  the  numerical  techniques  is  appropriate 
for  the  class  of  nonlinear  problems  associated  with  predicting  limit  cycle 
response  of  aerodynamic  surfaces  containing  discrete  structural  nonlinear¬ 
ities.  A  second  conclusion  of  this  study  is  that  the  more  refined  numerical 
simulations  yielded  improved  correlation  with  the  approximate  solution  of 
Reference  9,  when  compared  with  the  trapezoidal  integration  scheme. 

Based  on  computational  efficiency,  the  Runge-Kutta  integration  approach 
appears  most  attractive  for  these  type  problems.  The  normalized  cost  of  the 
three  techniques  (CPU  time,  connect  time,  etc.)  are  shown  in  Figure  16.  As  can 
be  seen,  the  Runge-Kutta  approach  is  the  most  efficient  technique.  The  Shanks 
technique  is  higher  order  than  the  Runge-Kutta  method  and  requires  signifi¬ 
cantly  more  calculations  per  time  step.  The  Adams-Moulton  procedure  is  a 
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predi ctor-corrector  scheme  and  spends  a  great  deal  of  time  "stopping  and 
starting. " 

In  spite  of  the  improved  correlation  obtained  with  the  refined  integration 
schemes,  there  remained  regions  where  comparisons  between  the  numerical 
simulation  results  and  the  approximate  solutions  were  inconsistent.  These 
inconsistencies  were  most  apparent  for  a  flexible  aerodynamic  control  surface 
containing  two  root  spring  preload  nonlinearities.  Additionally,  these 
inconsistencies  appeared  in  limited  regions  of  the  nonlinear  system  response. 
This  latter  point  is  illustrated  by  the  data  presented  in  Figure  17.  The  trend 
in  simulation  results  does  not  agree  with  the  second  order  approximate  solution 
in  the  region  of  dynamic  pressure  between  55  and  65  psi.  As  discussed  in 
earlier  paragraphs,  limited  analysis  of  the  simulation  results  indicate  that 
there  appears  to  be  a  change  in  the  nature  of  the  limit  cycle  response  in  this 
region. 

The  mechanism  of  the  nonlinear  system  response  described  in  the  preceeding 
paragraph  is  not  thoroughly  understood.  This  indicates  the  need  for  additional 
research  to  provide  an  understanding  of  this  nonlinear  response  phenomenon. 
Study  activities  are  needed  to  provide  explanations  for  the  prediction 
inconsistencies  noted  in  this  study.  Additional  numerical  simulations  should 
be  performed  for  regions  deemed  of  interest,  such  as  the  "flat  spot"  noted  in 
Figure  17.  The  influence  of  various  parameters  on  the  interpretation  of  these 
simulation  results  need  to  be  evaluated.  Significance  of  parameters  such  as 
simulation  run  time,  initial  condition  magnitude  and  combination,  and  RMS 
amplitude  determi nat. i on  need  additional  evaluation.  In  addition,  detailed 
analyses  of  the  frequency  content  of  the  numerically  predicted  system  response 
are  of  interest.  These  activities  will  develop  an  improved  understanding  of 
the  nonlinear  response  characteristics  of  an  aerodynamic  surface  containing 
discrete  structural  non  1 i near i t i es . 
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FIGURE  17  REGION  OF  CORRELATION  ANOMALY 
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APPENDIX  A 

AERODYNAMIC  SURFACE  CONFIGURATION 


Properties  of  the  Harpoon  missile  quick-attach  control  surface  were  used 
to  define  the  aerodynamic  surface  conf iguration  used  throughout  this  study. 
The  geometric  conf igurat ion  of  the  aerodynamic  surface  is  shown  in  Figure  Al. 
The  structural  nonlinearities  that  were  investigated  are  associated  with  root 
support  structure.  Presented  in  Figure  A2  are  the  inertia  properties  of  the 
aerodynamic  surface.  These  are  the  specific  terms  of  the  inertia  matrix  of 
Equation  (2).  The  first  two  rows  and  columns  of  the  inertia  matrix  are 
associated  with  rigid  root  roll  and  pitch  motions  while  the  last  two  diagonal 
elements  are  the  generalized  masses  of  the  aerodynamic  surface  modes.  The  off- 
diagonal  terms,  the  PF  quantities,  represent  the  inertia  coupling  between 
rigid  and  flexible  motions.  The  mode  shapes  associated  with  the  first  two 
aerodynamic  surface  cantilever  modes  are  given  in  Figure  A3.  These  modal  data 
were  used  when  investigating  a  flexible  surface  configuration. 
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(a)  Form  of  Inertia 

matrix 

0.1667  lb-sec^— in. 

PF01 

0.071  Ib-sec^-m. 

PP02 

0.058  lb-sec^-m. 

PP02 

4.220  x  10-4  lb-sec^/in. 

PF02 

4.295  x  10~ 4  Ib-sec^/in. 

(b)  Specific  inertia  terms 


-7.227  x  10~3 
-1.014  x  10~3 
-3.212  x  10“3 
2.134  x  IQ-3  |t 


FIGURE  A2  AERODYNAMIC  SURFACE  INERTIA  PROPERTIES 


APPENDIX  B 

SIMULATION  RESULTS  FOR  FLEXIBLE  AERODYNAMIC  SURFACE 


Presented  on  the  following  pages  are  plots  of  the  root  roll  and  root  pitch 
response  obtained  from  the  numerical  simulation  for  the  flexible  aerodynamic 
surface.  These  results  were  obtained  for  a  surface  having  two  preload 
non! ineari ties  with  the  following  characteristics: 

S<t>  =  S©  =  0.2  deg 
P<j>  =  Pe  =  0.1  deg 

In  addition,  the  uncoupled  root  roll  frequency  ( cl>q )  was  60  Hz  and  the  pitch 
frequency  (to*)  was  215  Hz  for  these  simulation  cases.  The  plots  on  the 
following  pages  formed  the  basis  for  obtaining  the  results  presented  in  Figures 
i2  and  13. 
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